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We adapt a variational procedure to calculate ground state properties of the Holstein model 
in the adiabatic limit. At strong coupling, this adaption leads to rapid convergence of results. 
The intermediate coupling regime is further handled with an adaptive algorithm. We also use semi- 
classically derived results for the adiabatic end-point, along with weak coupling perturbation theory. 
These establish weak and strong coupling (or large and small polaron, respectively) regimes in two 
dimensions or higher. As is well known, these are connected smoothly, but the cross-over becomes 
increasingly abrupt as the phonon frequency decreases. 

PACS numbers: 



I. INTRODUCTION 

There has been considerable work performed over the 
last two decades on the Holstein modeli. Interest in this 
model is fueled by the fact that it serves as the paradigm 
for electron-phonon interactions, much like the Hubbard 
model2, serves the same purpose for electron-electron in- 
teractions. While a considerable amount of this work has 
focussed on the many-electron problem, another subset 
has examined the single-electron, or polaron problem. A 
recent review of this work is available, for example, in 
Ref. 13 . 

In our opinion, the most promising numerical tech- 
nique for determining polaron properties in the thermo- 
dynamic limit is the variational procedure outlined by 
Trugman and coworkers^. With this method properties 
such as the ground state energy and the effective mass are 
readily obtained, in any dimension, over almost all pa- 
rameter regimes. One range of parameter space that has 
remained difficult, however, is near the adiabatic limit, 
which is what we address in this paper. The actual adia- 
batic limit was first treated by Kabanov and Mashtakov^; 
they found that in one dimension (ID), the electron re- 
tains polaronic character for all electron-phonon coupling 
strengths, while in two dimensions and higher there is a 
critical coupling strength, below which the electron be- 
haves in a 'free-electron-like' manner, and above which 
it is polaronic. At the same time, away from the adi- 
abatic limit the problem is known from numerical solu- 
tions to have a smooth crossover as a function of coupling 
strength (i.e. no abrupt transition), so it is of interest to 
pursue this crossover as the phonon frequency decreases 
towards zero. This was done to some degree in Refs. ii, 
but only for rather small lattices in one dimension. Our 
aim is to examine this limit using the Trugman varia- 
tional technique^ii^. 

The outline of the paper is as follows. In the next sec- 
tion we outline the model, and establish notation, etc. 
In Section HI we describe some refinements to the varia- 
tional method, and provide some illustrative examples to 
demonstrate the improvement in convergence. In Section 
IV we provide some numerical results as the adiabatic pa- 
rameter iOE/t approaches zero. Also provided are some 



perturbation theory results^, which can be reinterpreted 
to provide constraints for the numerical results. In Sec- 
tion V we show some results concerning the expected 
numbers of phonons in the ground state, which gives an- 
other indication of the difficulty of the adiabatic limit. 
Finally, we close with a summary. 



II. THE MODEL 

The model that most simply describes an electron in- 
teracting with optical phonons is the Holstein Hamilto- 
nian, given by 



i i 

which is also written 

+ WB^aJa,, (2) 



where cl (ci) creates (annihilates) an electron at site i 
(the spin label is suppressed) and = c|ci is the number 
density operator. The ion momentum pi, and displace- 
ment Xi are quantized via 



^ (aj - a,) , 



(3) 



where AI is the ion mass (we set Ti = 1) and a] (a^) creates 
(annihilates) a phonon at site i. The sum over i is over all 
sites in the lattice, whereas the sum over S is over near- 
est neighbors. Here, as the notation already suggests, 
we confine ourselves to nearest neighbor hopping only. 
The parameters are the hopping integral t, the phonon 
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frequency w^;, and the coupling of the electron to the os- 
cillator degrees of freedom, a. This parameter is the bare 
coupling between the electron and the ion; however, it is 
rarely used, and instead in the polaron literature the di- 
mensionless coupling constant q = — — ° is used. In 
the many-body literature, the dimensionless parameter 
A = '^'^^^ = j^ j-^^ ;y is used, where W = 2zt is the elec- 
tronic bandwidth for a cubic tight-binding model with 
coordination number z (z = 2,4,6 in 1,2,3 dimensions, 
respectively) . The parameter A has historical significance 
for the effective mass of degenerate electrons weakly cou- 
pled to phonons. Alternatively, and most useful in the 
strong coupling regime, the parameter g (or g^) was used 
in the Lang-Firsov transformatioiiii', and leads to a band 
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narrowing factor t ^ t* ^ te^^ in first order degenerate 
perturbation theory. We write all energy scales in terms 
of the hopping integral, t, which, hereafter is set to unity. 
So two dimensionless parameters which usefully charac- 
terize this problem'^ are w^/t, the adiabaticity parame- 
ter, and, A = 2g'^uJE /W . Actually, as we argue below in 
Section III, in two dimensions it is arguably more useful 
to use A = 2g2tj£;/(47rt) = g'^ uj e / {2iTt) , where l/(47rt) is 
the value of the non-interacting electron density of states 
at the bottom of the band (as opposed to the average 
density of states). 

As mentioned in the introduction this model has been 
most successfully analyzed using a refinement of the stan- 
dard Lanczos method due to TrugmanAtS-. Very accu- 
rate results can be obtained in any dimension-^" in al- 
most all parameter regime o'^'^'^" . A difficulty remains for 
moderately to strongly coupled systems with low adia- 
baticity parameter w^/i. For example, if one uses the 
Lang-Firsov transformatioi>ii to define the zeroth order 
strongly coupled wave function, then the average num- 
ber of phonons in the ground state can be readily deter- 
mined to be approximately g^ . For typical parameters in 
the moderately coupled regime (in one dimension), say 
loe = O.OSt, and A = 1.0, then g^ — 40, and this is 
the approximate number of phonons in the ground state. 
The Trugman procedure starts with a bare electron; on a 
moderate work station a feasible number of applications 
of the Hamiltonian is Nh = 22 (as remarked in Ref. 0), 
which produces a Hilbert space of order 10^. This process 
with Nh — 22 produces states that contain a maximum of 
22 phonons, and cannot possibly yield the correct ground 
state. 



III. REFINEMENT OF THE TRUGMAN 
METHOD 

We have examined two simple refinements to the Trug- 
man method^^ ; instead of starting with the bare electron 
state (properly extended throughout an infinite lattice), 
we first start with the state which is used as the unper- 



turbed state in the strong-coupling limit^iiii: 

\^) = e-s'/2^e'''=«'e-s'ilc||0), (4) 
i 

where the sum is over all lattice site o^'^'^"^ . As we shall see 
in what follows this speeds up convergence considerably 
in the strong coupling regime (either A>>loraj£;<< 
t). An example of the increased convergence is illustrated 
in Fig. [TJ Following Fig. 2 of Ref. [13 we show the 
fractional error Ajv = \{En — Em-i) / Em\ as a function 
of the number of states kept in the Hilbert space, for two 
sets of parameters, both of which have g^ = 20, using the 
bare electron as the starting state and using the strong- 
coupling solution as the starting state. There is a clear 
numerical advantage to using the latter. In Fig. [Ijc) 
we show the fractional error for a parameter regime near 
A ~ 1, where the strong coupling start is better even 
for values of A < 1. It is also clear that as A increases 
beyond the range of this figure, the refinement becomes 
increasingly useful. 

In pursuit of more severe disparate electron and 
phonon energy scales we found that even starting with 
the strong coupling solution resulted in slow convergence 
when A was of order unity. A remedy to this difficulty 
is the following procedure: start at large values of A, 
where convergence is readily obtained after a few itera- 
tions. Lower the value of A by a small amount, and use as 
a starting wave function the previous solution, truncated 
to include components with some minimal amplitude (so 
that a few hundred basis states at most are used) . Then 
converge the solution for this value of A, lower it, and 
continue the process until the desired range is covered. 
We have found this latter procedure to be the most ro- 
bust, particularly when the phonon frequency is much 
smaller than the electron hopping parameter, t. 

In Fig. [2]we show the ground state energy £'o vs. A for 
various values of the phonon frequency; this is in one di- 
mension. Fig. [3] shows similar results in two dimensions. 
It is clear that as the phonon frequency decreases, the 
crossover region near A ~ 1 (actually, the 'critical' value 
of A, only valid in the adiabatic limit, is closer to 0.55) 
becomes sharper. This is consistent with the result that, 
in the adiabatic limit, there is a transition from a small 
polaron state to a free electron-like state, in dimensions 
two and higher—. Nonetheless, as is known through other 
considerations^^, for any non-zero phonon frequency, the 
crossover is smooth. 

To summarize this section, we have obtained numeri- 
cally exact results for a wide range of parameters, by us- 
ing refinements to that used in Ref. [^Isi In particular, we 
obtain well converged results over all coupling strengths 
and for low phonon frequencies, we « t. The results 
for low frequencies in particular illustrate a rather abrupt 
crossover to a regime where multi-phonon processes are 
prevalent. To what extent they play a crucial role even 
at intermediate coupling strengths is the subject of the 
next section. 
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FIG. 1: (color online) Fractional error vs. number of states 
for (a) a very strong coupling case, and (b) a moderate cou- 
pling, but with itJEjt = 0.1. In both cases our refinement 
speeds up convergence considerably. In (c) we show the im- 
provement over a range of coupling strength near A = 1; for 
A > 1 (not shown) it is clear that the refinement leads to 
much better convergence. Note that the fractional error does 
not have to decrease monotonically with the number of states 
added. 



FIG. 2: (color online) Ground state energy i?o vs. A and 
m*/m vs. A for various phonon frequencies, in one dimension. 
There really is no special value of A singled out in these curves, 
consistent with the crossover phenomenon discussed in the 
text. 



IV. PERTURBATION THEORY 



Perturbation theory can be performed both from the 
weak and the strong coupling limits In that work we 
obtained, to second order (in 17), in one dimension, the 
self energy 



Eid(w -f iS) = 



Aa;£;sgn(a; — uje) 



(5) 
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FIG. 3: (color online) (a) Ground state energy Eq vs. A for 
various phonon frequencies, in two dimensions. In contrast to 
the ID results in Fig. [2] a 'special' value of A is now apparent 
— A « 0.55. However, for any non-zero phonon frequency 
the behaviour below and above this special value is smoothly 
connected. Only in the adiabatic limit does the behaviour 
change abruptly. 

(b) Expansion of the weak coupling regime showing the nu- 
merical results along side the perturbation theory results. 
Agreement is very good. 



FIG. 4: (color online) Effective mass m*/m vs. A for var- 
ious phonon frequencies, in two dimensions. Again, unlike 
the results in ID in Fig. [2] a 'special' value of A is clear — 
A f» 0.55. However, for any non-zero phonon frequency the 
behaviour below and above this special value is smoothly con- 
nected. Only in the adiabatic limit does the behaviour change 
abruptly. 

(b) Expansion of the weak coupling regime showing the 
numerical results alongside the perturbation theory results. 
Agreement is not as good as in Fig. 3. 



which leads to a ground state energy: 

This expression can be understood as follows: for very 
large frequency lue >> t there is a correction by a fac- 
tor 1 -|- A reminiscent of the mass renormalization for the 



electron mass in a Fermi liquid state. On the other hand, 
as the frequency becomes small the first order correction 
vanishes. In fact, however, the most significant effect of 
the phonon coupling to a single electron occurs for low 
phonon frequencies, while the effect disappears for high 
phonon frequency. This is most readily seen by examin- 
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ing the quasi-particle residue 



approximate form for the ground state energy is: 



^ ' 2t 



2\ UE (1 + ^)3/2- 



or the effective mass, defined as 



m/m* 



1 d'^E{k) I 



|fc=0- 



(7) 



(8) 



For a momentum independent self energy (as in the sec- 
ond order wealc coupling expansion) these are simply re- 
lated: m*/m = 1/ zq. The residue clearly approaches the 
non-interacting value, unity, as — >■ oo, while it di- 
verges as — )■ 0. This indicates a breakdown in (weak 
coupling) perturbation theory in this limit, which is con- 
sistent with the fact established in Ref. [zl that the elec- 
tron is polaron-like for all coupling strengths, i.e. there is 
an abrupt change in character aX g = Q. In fact, as estab- 
lished in Ref. [l[ for a two-site model, and in Ref. [7| , the 
effective mass diverges in the adiabatic limit for all cou- 
pling strengths (in ID), a limit which we now approach 
numerically in Fig. 2b. 

In two dimensions, as mentioned in Section II, we use 
A ~ bj E I i^T^t) ■ This actually uses the electron density 
of states at the bottom of the band, A^(— 4i) = l/(47r/:), 
instead of the average density of states that is commonly 
used, A'^avc(O) = l/(8i). The reason for this choice is 
that we are studying the one electron sector, so the most 
pertinent density of states is the one at the bottom. 

The self energy (in two dimensions (2D)) in weak cou- 
pling is given by 



I ^ ■X\ ^ 8iW-E 

S2d(w + 10) = K 



2 LO — LOe 



At 



uj — uje 



(9) 



where K{x) = J^^^ dO 



is the complete Elliptic 



^/l-ajsin^ e 

integral of the first kind. This leads to a ground state 
energy, which, in weak coupling, is: 
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We can take the derivative of Eq. © to obtain: 



A 



m* /m = 1 H 



1 



2 1 



-E 



1 



[1 + UJE I my 



(10) 



(11) 



LOEim 

where E{x) = /J^^^ dO^/ 1 — x sin^ is the complete El 
liptic integral of the second kind. We have used, ^"^^^-^ = 

K{x 



2x 



E{x) 
l-x 



More familiar expressions are available, for cases when 
the arguments of the complete elliptic integrals are 
close to unity. This occurs for << t. Using 
\\mm^iK{m) — iln(16/mi)yi^ where mi = 1 — m, an 



-4t(l+^"^ 



1 



4 t l + LOE/m 
From Eq. pT|) we obtain: 



In 



2t \+ujE/m 



(12)" 



m* /m = 1 



2 16t \ijJE J \ \ t J \u}E 



(13) 

Note that as uJE/t — > 0, the ground state energy ap- 
proaches the non-interacting value, while the result for 
the effective mass approaches the one derived in the con- 
tinuum limit by Cappelluti et alJ^, and the mass en- 
hancement is half that expected when Ep is large. 

In Fig. 3a, we show the ground state energy of the 2D 
Holstein model as a function of A for a variety of phonon 
frequencies; we also show the result in the adiabatic limit 
as ws/i — 7- 0. For the latter case, we adopted the iter- 
ative method described in Ref. [7,^, and used Lanczos 
diagonalization for the electronic portion. The abrupt 
transition occurs because we do not assume Bloch's the- 
orem, and translational invariance is broken for suffi- 
ciently strong coupling. For non-zero phonon frequency 
we note the trend that as — > 0, the crossover from 
free-electron-like behaviour to polaronic behaviour be- 
comes more abrupt, though it is always smooth. 

In Fig. 3b we show an expanded region in the weak 
coupling regime, where the perturbation theory results 
are also plotted. Note that they are quite accurate for 
all frequencies shown. 

In Fig. 4a, we show the electron effective mass for the 
same parameters as in Fig. 3. In strong coupling the 
effective mass grows rapidly with coupling strength, as 
shown. However, this increase is even more pronounced 
as the phonon frequency decreases, until, as the adiabatic 
limit is approached, the increase becomes very nearly 
abrupt above a 'critical' coupling strength, as determined 
through the adiabatic calculation. At weak coupling, 
the effective mass is unity for large u)e (not shown). 
As the phonon frequency decreases, the effective mass 
grows; however, for smaller phonon frequencies the ef- 
fective mass will decrease again as the phonon frequency 
decreases (as can be seen from the cases shown). Both of 
these trends conspire to make the crossover more abrupt 
as the phonon frequency approaches zero. 

In Fig. 4b we show an expanded region in the weak 
coupling regime (no log scale), where the perturbation 
theory results are also plotted. The results are certainly 
not as accurate as the ground state energy; however, the 
inversion with phonon frequency noted above is clearly 
obtained. 
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FIG. 5: (color online) (a) Mean plionon number vs. cou- 
pling strength for various plionon frequencies (in 2D). Note 
the increasingly abrupt behaviour as the phonon frequency 
decreases. 

(b) Mean phonon number as a function of wave vector. Even 
for very small coupling strength there is an abrupt increase 
when the phonon frequency is small enough. Explanation is 
provided in the text, and is confirmed by (c) where the energy 
as a function of wave vector is plotted for the same parameters 
as in (b). 



V. MEAN PHONON NUMBER AND 
DISPERSION ANOMALIES 

Finally, we briefly examine the expectation of the num- 
ber of phonons in the ground state, and the impact on 
the electronic dispersion relation. Restricting ourselves 
to two dimensions, we plot, in Fig. 5a, the mean phonon 
number vs. coupling strength in the intermediate cou- 
pling regime for several phonon frequencies. The same 
trend as seen in Fig. 4 is apparent — beyond a 'special' 
coupling strength the mean phonon number grows very 
abruptly from near zero to some value, N^, after which 
it continues to grow gradually as the coupling strength 
increases. The actual value of Ng is close to the central 
value of the Poisson distribution as predicted by strong 
coupling perturbation theory^. 

In Fig. 5b, we shows numerical results of the mean 
phonon number as a function of total momentum k^, of 
the electron-phonon system, for a few values of phonon 
frequency, and for a very low value of coupling (so that 
the results are well converged). Despite this small value 
of coupling, convergence is difficult because we use uje 
as small as O.Olt. We apply our self-adaptive Lanczos 
method by first converging the results for some high mo- 
mentum (say, fca;=0.3 — we keep ky = 0), and then lower 
the value of kx in small increments, and converge the cal- 
culation at each step, until we finally reach the desired 
end-point {k^ — 0)- As Fig. 5b illustrates, for sufficiently 
small phonon frequency, the mean phonon number shows 
a sharp increase from close to zero to nearly unity at some 
wave- vector, say kc- The reason for this is that the energy 
difference with the ground state will eventually exceed a 
value of order w^; at this point it becomes energetically 
more favourable to use the zero momentum state (with 
much lower energy), and simply excite a phonon with 
the required momentum. Confirmation of this explana- 
tion is provided in Fig. 5(c), where the dispersion flat- 
tens abruptly beyond kc, when the energy exceeds that 
of the ground state by an amount approximately equal 
to coe- It retains this value because phonon momenta of 
any value are available with the same energy. 



VI. SUMMARY 

We have implemented an adaptation to the variational 
method first suggested by Trugman, specifically to han- 
dle the adiabatic regime. In strong coupling our starting 
point leads to immediate convergence, while in the in- 
termediate coupling regime a 'stepping-down' procedure 
allows for good convergence. Even in weak coupling, if 
the phonon frequency is significantly lower than the hop- 
ping parameter, our adaptive method is helpful, if not 
necessary. 

By determining ground state properties as a function 
of decreasing phonon frequency we have established a 
connection between numerical results at small but non- 
zero phonon frequency, and adiabatic limit results ob- 
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tained by using a semi-classical iterative procedure. It is 
clear that in one dimension no weak coupling perturba- 
tion regime exists, while in two dimensions (and higher) 
a definite weak coupling regime exists, and results de- 
rived within perturbation theory agree well with numer- 
ical results down to very low frequencies. Finally, as the 
phonon frequency decreases, more and more phonons are 
present in the ground state wave function, and these lead 
to anomalies in the electron dispersion relation. 

It is somewhat ironic that the most significant pola- 
ronic effects occur in the adiabatic regime, as uJE/t — > 0. 
This is where weak coupling perturbation theory breaks 
down completely. The second order result, which is sim- 
ply the so-called non-crossing approximation, fails to cap- 
ture the rapid onset of multi-phonon excitations that 
form an integral part of the ground state wave func- 
tion, as exemplified, for example, in Eq. Q; this is a 
breakdown that, for example, AlexandroviS has repeat- 
edly emphasized. At the same time, the so-called Migdal 
approximatior>i^, so key to the Eliashberg theory of su- 
perconductivity, is valid only in this limit. One then re- 
quires an understanding of how polaronic effects become 
minimized as more and more electrons are included in 



the problem. Apparently Pauli blocking plays an impor- 
tant role in mitigating the multi-phonon processes that 
constitute a single polaron. Future work will attempt to 
investigate this cross-over from polaron to weak coupling 
behaviour. 



Acknowledgments 

This work was supported in part by the Natural 
Sciences and Engineering Research Council of Canada 
(NSERC), by ICORE (Alberta), by Alberta Ingenuity, 
and by the Canadian Institute for Advanced Research 
(ClfAR). We thank Stuart Trugman for helpful corre- 
spondence in the early part of this work. FM is grateful 
to the Aspen Center for Physics, where some of this work 
was done (summer, 2007). 

^ present address: Dept. of Oncology, University of 
Alberta, Edmonton, AB, Canada T6G 1Z2 
^ present address: Dept. of Mathematics, University of 
British Columbia, Vancouver, BC, Canada 



1 T. Holstein, Ann. Phys. (New York) 8, 325 (1959). 
^ J. Hubbard, Proc. Roy. Soc. A 276, 238 (1963); ibid., 277, 
237 (1964). 

^ H. Fehske and S.A. Trugman, in Polarons in Advanced 
Materials edited by A. S. Alexandrov, Springer Series in 
Material Sciences 103 pp. 393-461, Springer Verlag, Dor- 
drecht (2007). 

* A.S. Alexandrov, in Polarons in Advanced Materials, 
edited by A.S. Alexandrov, Springer Series in Materials 
Science, 103, pp. 257-310, Springer Verlag, Dordrecht 
(2007). 

^ S.A. Trugman, in Applications of Statistical and Field The- 
ory Methods to Condensed Matter, edited by D. Baeriswyl, 
A.R. Bishop, and J. Carmelo (Plenum Press, New York, 
1990). 

® J. Bonca, S.A. Trugman, and I. Batisti'c, Phys. Rev. B60, 
1633 (1999). 

V.V. Kabanov and O.Yu Mashtakov, Phys. Rev. B47, 6060 
(1993). 

* A.S. Alexandrov, V.V. Kabanov, and D.K. Ray, Phys. Rev. 
B49, 9915 (1994). 

^ F. Marsiglio, Physica C244 21, (1995). 
° L-C. Ku, S.A. Trugman, and J. Bonca, Phys. Reb. B65, 
174306 (2002). 



" I.G. Lang and Yu. A. Firsov, Sov. Phys. JETP16, 1301 
(1963); Sov. Phys. Solid State 5 2049 (1964). 
A preliminary report of this work was given in Z. Li, D. 
Baillie, C. Blois, and F. Marsiglio, Bull. Am. Phys. Soc. 
54, H13.9 (2009). 

In the implementation of Trugman's method the use of a 
Bloch wave is implicit; nonetheless we include it explicitly 
in Eq. (|4]). Furthermore, in this representation we utilize 
phonon number states to reconstruct the exponential up 
to some high cutoff. 

Two other possibilities for improving the original Trugman 
method, which we learned about after this work was com- 
pleted, can be found in J. Bonca et al., Phys. Rev. B77, 
054519 (2008), and A. Alvermann et al. arXivilOOl. 24821 
H. Lowen, Phys. Rev. B37, 8661 (1988). 
M. Abramowitz and LA. Stegun, Handbook of Mathemat- 
ical Functions (Dover, New York, 1972). 
E. Cappelluti, C. Grimaldi, and F. Marsiglio, Phys. Rev. 
B76, 085334 (2007). 

A.B.Migdal, Zh. Eksp. Teor.Fiz. 34, 1438 (1958) [Sov. 
Phys. JETP 7, 996 (1958)]. 

A.S. Alexandrov, Europhys. Lett. 56, 92 (2001). 



